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Experimentally it is known that the bulk modulus, K, and shear modulus, ^, of a granular assembly 
of elastic spheres increase with pressure, p, faster than the p^^^ law predicted by effective medium 
theory (EMT) based on Hertz-Mindlin contact forces. Further, the ratio K/ ^ is found to be roughly 
pressure independent but the measured values are considerably larger than the EMT predictions. To 
understand the origin of these discrepancies, we have undertaken numerical simulations of a granular 
assembly of spherical elastic grains. Our results for K{p) and ^{p) are in good agreement with the 
existing experimental data. We show, also, that EMT can describe their pressure dependence if 
one takes into account the fact that the number of grain-grain contacts increases with p. Most 
important, the ajfine assumption (which underlies EMT), is found to be valid for K{p) but to 
breakdown seriously for /u(p). This explains why the experimental and numerical values of ii(p) are 
much smaller than the EMT predictions. 
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The study of sound propagation and nonlinear elastic- 
ity in unconsolidated granular matter is a topic of great 
current interest In the simplest experiments, a pack- 
ing of glass beads is confined under hydrostatic conditions 
and the compressional and shear sound speeds, Vp and 
Va, are measured as functions of pressure, p 1,1 . In the 
long- wavelength limit, the sound speeds are related to the 
elastic constants of the aggregate: Vp = \J{K + A/i^)/ p* 
and Ws = \/ isl p* , where p* is the system's density. In 
acoustic measurements were made on 



a recent Letter [E 
bead packs under uniaxial stress and it was suggested 
that long wavelength compresional waves can be de- 
scribed in terms of an effective medium. Thus, it would 
of great value to have a reliable EMT to describe sound 
propagation as a function of applied stress. However, our 
analysis, together with the work of others, raises serious 
question about the validity of the generally accepted the- 
oretical formulation. The EMT predicts that K and 
/i both vary as p^^^, and that the ratio K/ p is a constant 
(independent of pressure and coordination number) de- 
pendent only on the Poisson's ratio of the material from 
which the individual grains are made. 

Experimentally (see Fig. m, the bulk and shear moduli 
increase more rapidly than p^^'^ and the values of K / p are 
considerably larger than the EMT prediction. These dis- 
crepancies between theory and experiment could be due 
to the breakdown of the Hertz-Mindlin force law at each 
grain contact as proposed in ||^ for the case of metallic 
beads with an oxide layer, and in Q for grains with sharp 
angularities. Alternatively, they could be associated with 
the breakdown of some of the assumptions underlying 
the EMT, for example, that the number of contacts per 
grain is pressure independent, which may not be the case 
as several authors have suggested 0|8|. 

In this Letter we report calculations of K(p) and 
p{p) based on granular dynamics (GD) simulations using 



the Discrete Element Method developed by Cundall and 
Strack Here it is assumed at the outset that one 

has an assembly of spherical soft grains which interact via 
the Hertz-Mindlin force laws. We find good agreement 
with the existing experimental data, thus confirming the 
validity of Hertz-Mindlin contact theory to glass bead ag- 
gregates. Further, we can explain the two problems with 
EMT described above. First, if the calculated increase of 
the average coordination number with p is taken into ac- 
count, the modified EMT gives an accurate description of 
the bulk modulus found in the simulations, K{p)] for p{p) 
we obtain a curve whose shape is in good agreement with 
the simulation data but whose values are seriously offset 
therefrom. Second, the EMT makes the affine assump- 
tion in which the motion of each grain is specified simply 
in terms of the externally applied strain (see below). We 
show that while the affine assumption is approximately 
valid for the bulk modulus, it is seriously in error for the 
shear modulus; this is why the EMT prediction of K/ p 
differs significantly from the experimental value. 

Numerical Simulations. At the microscopic level the 
grains interact with one another via (1) non- linear Hertz 
normal forces and (2) friction generated transverse forces. 
The normal force, /„, has the typical 3/2 power law de- 
pendence on the overlap between two spheres in contact, 
while the transverse force, ft depends on both the shear 
and normal displacements between the spheres ||l^ . For 
two spherical grains with radii i?i and i?2: 



(la) 



(lb) 



Here R = 2i?ii?2/(i?i 

W = (l/2)[(i?i +i?2) - 



+ i?2), the normal overlap is 
2^1 ^2^21] > 0, and xi, X2 are 
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the positions of the grain centers. The normal force acts 
only in compression, fn = when w < 0. The vari- 
able s is defined such that the relative shear displace- 
ment between the two grain centers is 2s. The prefactors 
C„ = 4G/(1 - ly) and Ct = 8G/{2 ~ v) are defined in 
terms of the shear modulus G and the Poisson's ratio v 
of the material from which the grains are made. In our 
simulations we set G = 29 GPa and v = 0.2. We assume 
a distribution of grain radii in which i?i = 0.105 mm for 
half the grains and i?2 = 0.095 mm for the other half. 
Our results are, in fact, quite insensitive to the choice of 
distribution. We also include a viscous damping term to 
allow the system to relax toward static equilibrium p2[ . 
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FIG. 1. Pressure dependence of the elastic moduli from GD 
(•), experiments [Domenico (□) Yin (O) Q, and our ex- 
periments (O)]; ^-iid theory (dashed line): (a) bulk modulus, 
and (b) shear modulus. 

Our calculations begin with a numerical protocol de- 
signed to mimic the experimental procedure used to pre- 
pare dense packed granular materials. In the experiments 
the initial bead pack is subjected to mechanical tapping 
and ultrasonic vibration in order to increase the solid 
phase volume fraction, (ps- The simulations begin with 
a gas of 10000 spherical particles located at random po- 
sitions in a periodically repeated cubic unit cell approx- 
imately 4 mm on a side. At the outset, the transverse 
force between the grains is turned off (Ct = 0). The sys- 
tem is them compressed slowly until a specified value of 
ips is attained (see dashed lines in Fig. ^). The compres- 
sion is then stopped and the grains are allowed to relax. 
If the compression is stopped before reaching the criti- 
cal volume fraction, 4>s ~ 0.64, corresponding to random 
close packing (RCP), the system will relax to zero pres- 
sure and zero coordination number, since the system can- 
not equilibrate below RCP. The compression is then con- 
tinued to a point above the critical packing fraction and 
the target pressure is maintained with a "servo" mecha- 
nism 1^ which constantly adjusts the applied strain until 
the system reaches equilibrium. Because there are no 
transverse forces, the grains slide without resistance dur- 
ing the relaxation process and the system reaches the 
high volume fractions found experimentally. 

The simulated granular aggregate relaxes to equilib- 
rium states in which the average coordination number, 
< Z{p) >, increases with pressure as seen in Fig. ||. 



For low pressures compared with the shear modulus of 
the beads < Z >w 6, while in two dimensions the same 
preparation protocol gives < Z >« 4. Such low coor- 
dination numbers can be understood in terms of a con- 
straint argument for frictionless rigid balls |]l3|,^ , which 
gives < Z 2D, where D is the dimension. These val- 
ues should be valid in the limit of low pressure when the 
beads are minimally connected near RCP |lj] (or in the 
rigid ball limit G — > oo). For large values of the confining 
pressure more grains are brought into contact, and the 
coordination number increases . Empirically, we find 
that 



< Z{p) >= 6 
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Comparison with Experiment. Consider, now, the cal- 
culation of the elastic moduli of the system as function 
of pressure. Beginning with the equilibrium state de- 
scribe above, we first restore the transverse component 
of the contact force interaction. We then apply a small 
perturbation to the system and measure the resulting re- 
sponse. The shear modulus is calculated in two ways, 
from a pure shear test, fi = (l/2)Acri2/Aei2, and also 
from a biaxial test, fj, = (A(T22 — Ao'ii)/2(Ae22 — Aen). 
The bulk modulus is obtain from a uniaxial compression 
test, K + 4/3/i = AcTii/Aeii. Here the stress, Cy, is 
determined from the measured forces on the grains 
and the strain, e^j, is determined from the imposed di- 
mensions of the unit cell. 
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FIG. 2. Coordination number versus pressure. The 
dashed line is a fit according to Eq. ^. 

In Fig. [l] our calculated values of the elastic moduli as 
a function of pressure are compared with EMT and with 
experimental data. Because there is a considerable de- 
gree of scatter in the experimental results we performed 
our own experiments with standard sound propagation 
techniques. A set of high quality glass beads of diam- 
eter 45/im are deposited in a flexible container of 3 cm 
height and 2.5 cm radius. Transducers and a pair of lin- 
ear variable differential transformers (for measurement of 
displacement) were placed at the top and bottom of the 
flexible membrane, and the entire system was put into a 
pressure vessel filled with oil. Before starting the mea- 
surements, a series of tapping and ultrasonic vibrations 
were applied to the container in order to let the grains 
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settle into the best possible packing. We then applied 
confining pressures ranging from 5 MPa to 100 MPa. The 
pressure was cycled up and down several times until the 
system exhibited minimal hysteresis. At this point shear 
and compressional waves were propagated by applying 
pulses with center frequencies of 500 KHz. The sound 
speeds and corresponding moduli were obtained by mea- 
suring the arrival time for the two sound waves. 

From Fig. |l] we see that our experimental and nu- 
merical results are in reasonably good agreement. Also 
shown are measured data from Donienico Q and Yin . 
Clearly, the experimental data are somewhat scattered. 
This scatter reflects the difficulty of the measurements, 
especially at the lowest pressures where there is signifi- 
cant signal loss. Nevertheless, our calculated results pass 
through the collection of available data. Also shown in 
Fig. Ilfare the EMT predictions H 
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The EMT curves are obtained using the same parameters 
as in the simulations; we also set Z = 6 and 4>s = 0.64, 
independent of pressure. At low pressures we see that K 
is well described by EMT. At larger pressures, however, 
the experimental and numerical values of K grow faster 
than the p^/^ law predicted by EMT. The situation with 
the shear modulus is even less satisfactory. EMT overes- 
timates fj,(j)) at low pressures but, again, underestimates 
the increase in /i(p) with pressure. 
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FIG. 3. Shear modulus and corrected EMT taking into ac- 
count the pressure dependence of < Z{p) > from Fig. |^ as 
well as ^s(p). 

To investigate the failure of EMT in predicting the cor- 
rect pressure dependence of the moduli, we plot in Fig. ^ 
the shear modulus divided by p^^"^. For such a plot, EMT 
predicts a horizontal straight line but we see that the nu- 
merical and experimental results are clearly increasing 
with p. We can explain this behavior by modifying Eq. 
(3b) to take account of the pressure dependence of the 
coordination number < Z{p) > from Fig. |^ and the 
pressure dependence of (psip) (which is a much smaller 



effect) . The modified EMT is also plotted in Fig. ^ and 
we see that it predicts the same trend with pressure as 
the simulations. The experimental data also seem to be 
following this trend but more data over a larger pressure 
range are clearly needed. Not shown in Fig. |^ is a sim- 
ilar analysis of K{p) but the result is that the modified 
EMT is in essentially exact agreement with our numer- 
ical simulations. It is for this reason that we focus on 

Another way of seeing the breakdown of EMT is to 
focus on the ratio K/ii. According to Eqs. (||), K/fi ~ 
5/3(2 — i')/{5 — 4i/), independent of pressure, a value 
which depends only on the Poisson's ratio of the bead ma- 
terial. The experiments give K/ fj, « 1.1 — 1.3. Our simu- 
lations give K/^ w 1.05±0.1 in good agreement with ex- 
periments. EMT predicts K/ii = 0.71, if we take = 0.2 
for the Poisson's ratio of glass. [The EMT prediction is 
quite insensitive to variations of v; K/fi — 0.71 ±0.04 for 
ly = 0.2 ±0.1.] 

To understand why /i is over predicted by EMT we 
must examine the role of transverse forces and rotations 
in the relaxation of the grains. [These effects do not play 
any role in the calculation of the bulk modulus.] Sup- 
pose we re-define the transverse force by introducing a 
multiphcative coefficient a, viz: Aft — aCt^RwY^'^As; 
with a = I we recover our previous results. To quantify 
the role of the transverse force on the elastic moduli, we 
calculate K{a) and /i(a) at a given pressure p = 100 KPa 
[Fig. I^a]. This pressure is low enough that the chang- 
ing number of contacts is not an issue. Surprisingly, the 
shear modulus becomes negligibly small as a — > 0. As 
expected, K is independent of the strength of the trans- 
verse force. To compare with the theory we also plot the 
prediction of the EMT, Eq. (^ in which is rescaled by 
aCf. We see that the EMT fails in taking into account 
the vanishing of the shear modulus as a — > 0. However it 
accurately predicts the value of the bulk modulus, which 
is independent of a. 

There are two main approximations in the EMT: (1) 
All the spheres are statistically the same, and it is as- 
sumed that there is an isotropic distribution of contacts 
around a given sphere. (2) An afhne approximation is 
used, i.e., the spheres at position Xj are moved a dis- 
tance Sui in a time interval 6t according to the macro- 
scopic strain rate Sij by 6ui = iijXjSt. The grains are 
always at equilibrium due to the assumption of isotropic 
distribution of contacts and further relaxation is not re- 
quired. 

In the GD calculation of the shear modulus an affine 
perturbation is first applied to the system. The shear 
stress increases (from A to B in Fig. ^o) and the grains 
are far from equilibrium since the system is disordered. 
The grains then relax towards equilibrium (from B to C), 
and we measure the resulting change in stress from which 
the modulus is calculated. To better understand the ap- 
proximations involved in the EMT, suppose we repeat 
the GD calculations taking into account only the affine 
motion of the grains and ignoring the subsequent relax- 
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ation. The resulting values of the moduli are plotted in 
Fig. ^a as open symbols and we see that the moduli cal- 
culated this way are very close to the EMT predictions. 
Thus, the difference between the GD and EMT results 
for the shear modulus lies in the non-afhne relaxation of 
the grains; this difference being largest when there is no 
transverse force. By contrast, grain relaxation after an 
applied isotropic affine perturbation is not particularly 
significant and the EMT predictions for the bulk modu- 
lus are quite accurate. 
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FIG. 4. (a) K{a), and /i(a) versus a for a fixed p — 100 
KPa. (b) Relaxation of the shear stress (B^C) after an affine 
motion (A^B) in the calculation of the shear modulus. 

The surprisingly small values of fi found as a ^ 
can be understood as a melting of the system that oc- 
curs when the system is close to the RCP fraction. This 
fluid like behavior (when Ct — 0) is closely related to the 
melting transition seen in compressed emulsions [T^ ] and 
foams . At the RCP fraction the system behaves like a 
fluid with no resistance to shear. By contrast, molecular 
dynamics simulations of glasses, in which the atoms inter- 
act by purely longitudinal forces, predict non-vanishing 
shear speeds [|l8|. The crucial difference between these 
two systems is the local coordination of the particles. In 
the granular system, the coordination number near RCP 
(where the balls are only weakly deformed) is < Z >= 6; 
the system is at its minimal coordination number. In 
glasses, however, the number of neighbors is closer to 10 
and the motion of the grain is highly constrained. 

Conclusions. Our GD simulations are in good agree- 
ment with the available experimental data on the pres- 
sure dependence of the elastic moduli of granular pack- 
ings. They also serve to clarify the deficiencies of EMT. 
Grain relaxation after an infinitesimal afhne strain trans- 
formation is an essential component of the shear (but 
not the bulk) modulus. This relaxation is not taken into 
account in the EMT. In the limit a a packing of 
nearly rigid particles responds to an external isotropic 



load with an elastic deformation and a finite K. By con- 
trast, such a system cannot support a shear load (/z — > 0) 
without severe particle rearrangements. This may in- 
dicate a "fragile" state of the system [|l^ where inter 
particle forces are organized along "force chains" (stress 
paths carrying most of the forces in the system) oriented 
along the principal stress axes. Such fragile networks 
support, elastically, only perturbations compatible with 
the structure of force chains and deform plastically oth- 
erwise. Clearly, there is a need for an improved EMT; re- 
cent work on stress transmission in minimally connected 
networks |]l3|,^ may provide an alternative formulation 
and allow to describe properly the response of granular 
materials to perturbations. 
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